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The QCD evolution of both unpolarized and polarized gen- 
eralized parton distributions (GPDs) to next-to-leading order 
(NLO) accuracy is presented, in both the DGLAP and ERBL 
regions, for two appropriately symmetrized input distribu- 
tions based on conventional parton density functions. To il- 
lustrate the relative size of the NLO corrections a comparison 
is made with leading order evolution of the same distribu- 
tions. For the first time, NLO results are given for both small 
and large values of the skewedness parameter, £ = xtj, i.e. for 
all of the kinematic range relevant to HERA and HERMES. 



I. INTRODUCTION 

Generalized Partons Distributions (GPDs) [1-4] are 
generic two-parton correlation functions of nucleons that 
provide the boundary conditions for the calculation of 
various hard, exclusive, diffractive processes in QCD. For 
such processes, which are characterised by a suitable hard 
scale (for example a large photon virtuality, Q 2 , in Deep 
Inelastic Scattering (DIS) experiments), the incoming nu- 
cleon remains intact and at high energies is well separated 
in rapidity from the rest of the final state, i. e. the diffrac- 
tively produced low mass particle (e.g. a meson or real 
photon). A perturbative QCD analysis is possible pro- 
vided that a factorization theorem is proved which illus- 
trates the separation of short and long distance physics. 
The QCD amplitude is then a sum of terms involving a 
convolution of a hard scattering coefficient (calculated to 
a given order in a s ) with a GPD (specified to a related 
accuracy in logarithms) and possibly a second convolu- 
tion with a distribution amplitude, in the case of the 
production of a hadronic final state. 

The universal, non-perturbative GPDs obey renormal- 
ization group equations (RGE) and are formally defined, 
up to power suppressed contributions, by Fourier trans- 
forms of non-local, renormalizcd, light-cone operators 
sandwiched between nucleon states of unequal momen- 
tum (the corresponding distributions in inclusive DIS in- 
volve equal momenta for the incoming and outgoing nu- 
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cleon). In common with conventional DGLAP parton 
distribution functions, GPDs cannot be calculated from 
first principles at present due to the confinement prob- 
lem, but they can be pinned down through a global RGE 
analysis of all available experimental data. The simplest 
example of such a hard, exclusive, diffractive process is 
deeply virtual Compton scattering (DVCS) [1,2,5,6], for 
which there is recent first data available from HERA [7,8] , 
HERMES [9] and JLAB [10]. 

The RGE for the GPDs involve kernels which have 
been computed up to next-to-leading order (NLO) in per- 
turbation theory so far [11]. The solutions of these NLO 
evolution equations are necessary to obtain accurate pre- 
dictions for various experimentally accessible processes. 
In fact the DVCS amplitude itself can be directly ac- 
cessed experimentally through the interference of DVCS 
with the Bethe-Heitler process [7-9]. In this paper we 
present numerical solutions to the NLO QCD evolution 
of the GPDs for large (HERMES) and small (HERA) 
skewedness which corresponds to large and small Bjorkcn 
x in the particular representation used. As input for the 
NLO evolution, we use a model based on double distribu- 
tions (DDs) [12,13] incorporating for comparison two re- 
lated sets of conventional NLO unpolarized and polarized 
parton distribution functions (PDFs), MRSA [14] and 
Gehrmann-Stirling ('GS(A)') [15], and GRV98 [16] and 
GRSV00 [17], respectively. The resultant input GPDs, 
by design, respect all the known symmetries and proper- 
ties of the GPDs, which are preserved under evolution. 
This latter condition constitutes a stringent test of any 
proposed solutions, which is passed by our numerical so- 
lutions to very good accuracy. 

This paper presents the first results of the numerical 
solution to the RGEs at NLO accuracy [18-20] for any x 
(see, e.g., [21,22] for leading order skewed evolution) and 
is organised as follows. In section II we give the formal 
definitions of the GPDs and discuss their symmetry and 
spectral properties. Section III describes our choice of in- 
put distributions to the evolution equations outlined in 
section IV. We present our numerical results, including 
a comparison between leading order and next-to-leading 
order skewed evolution and a comparison with conven- 
tional next-to-leading order evolution, in section V, and 
conclude in section VI. Our evolution code will be made 
available to the community in the near future via the 
internet [23]. 
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II. DEFINITIONS, SYMMETRIES AND 
SPECTRAL REPRESENTATIONS OF GPDS 

Amplitudes for quark and gluon correlators of unequal 
momentum nucleon states may be defined in a number of 
ways. We choose a definition which treats the initial and 
final state nucleon momentum {p,p' , respectively) sym- 
metrically by involving parton light-cone fractions with 
respect to the momentum transfer, r = p — p', and the 
average momentum, P = (p + p')/2. The flavour singlet 
and non-singlet (S,NS) quark, and the gluon (G) matrix 
elements of the non-local operators, involving a light-cone 
vector (z 2 = 0), are defined by 

2M^ s (P-z,r-z,t) = 
(N{P_)\i> a (- Z -)Vz^ a (^) \N(P + )) + 

(N{P-)^ a (^)Vzi> a (- Z -) \N{P + )), 
2M S (P -z,r-z,t) = 

J2(N (P-) \i> a (-f ) Vz^ a (|) |JV (P+)> - 

a 

{N{P-)\i> a (l)Vz^ a (- Z -) \N(P + )), 
2M G (P ■ z,r- z,t) = 

z»z»{N (P_) \G, P (-|) VG1 (|) \N (P + )) + 
z»z v {N (P_) \G PP (|) VG'l (-|) \N (P + )}, 

(1) 

where P+._ = P±§, t = r 2 is the four-momentum trans- 
fer, z — 7 M 2 M , P • z and r ■ z are dimensionless Lorentz 
scalars, a = u,d, s, c... is a flavour index and the symbol 
V represents the usual path ordered exponential. The 
spin and isospin dependence have been suppressed for 
convenience. Polarized matrix elements are defined in 
a similar fashion, but with z — > z-y 5 for the quark and 
with the implicit metric g^ v — > ie? v v for the gluon. 
The unpolarizcd matrix elements, and thus their associ- 
ated GPDs, have definite symmetry properties. G-parity 
gives (suppressing the i-dependence) 



M NS (P ■ z,r ■ z) 



M S ' G {P -z,r-z) = M s ' u (-P -z,-r- z) 



-M NS (-P ■ z,-r- z) 

S,Gf 



(2) 



The opposite properties are observed for the polarized 
case. Hermitian conjugation gives for all polarized and 
unpolarized species (i = S, NS, G) 



M\P ■z,r-z) = M i (P-z,-r-z). 



(3) 



The matrix elements can be most generally represented 
by a double spectral representation with respect to P ■ z 



and r • z [1 ,3,24] as follows: 



M 



NS 



r l _ r l-\x\ 

/ dxe~ ixPz / dye 

J-l J-l+\x\ 



-iyr-z/2 , 



U'zUF NS {x,y,t) + 



K NS (x,y,t) 



M G — 



U'zU 



2rriN 

dxe~ lxP - z / dye- iyr - z ' 2 

-1 J-l + \x\ 

P-Z r , 



+ 4^ K ^ 

+ U'Ur ■ z j dye- tyr - z ' 2 D G {y, t) , 



M 



dxe 



-ixP-z 



l-|x| 



1 + M 



dye 



-iyr-z/2 „ 



U'zUF s (x,y,t) + 



2m n 



K s (x,y,t) 



+ 



U'Ur-z ^ dye- iyr - z/2 D{y,t), 



(4) 



where U and U are nucleon spinors. Note that in ac- 
cordance with the associated Lorentz structures, the P's 
correspond to helicity non-flip and the K's to helicity flip 
amplitudes, and are collectively known as double distri- 
butions. Henceforth, for brevity, we shall only discuss 
the helicity non-flip piece explicitly. However, the helic- 
ity flip case is exactly analogous. The D-terms in the last 
two lines permit non-zero values for the singlet M s and 
M G in the limit P ■ z — > and r-z^O, which is allowed 
by their evenness in P • z (cf. eq.(2)). Conversely, M NS 
is required by its oddness under P ■ z, for any r • z, to be 
zero [24]. 

By making a particular choice of the light-cone vec- 
tor, z* 1 , as a light-ray vector (so that in light-cone vari- 
ables, z± = zo±z 3 , only its minus component is non-zero 
z^ = (0, z_,0)) one may reduce the double spectral rep- 
resentation of eq.(4), defined on the entire light-cone, to 
a one dimensional spectral representation, defined along 
a light ray, depending on the skewedness parameter, £, 
defined by 



£ = r-z/2P-z = r + /2P + . 



(5) 



The resultant GPDs are the off-forward distribution func- 
tions (OFPDFs) introduced in [1,2]: 



l-|x'| 



dy'8(x' + Zy' -v)F(x',y',t), (6) 



where v G [—1,1]. In terms of individual flavour decom- 
position the singlet, non-singlet and gluon distributions 
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are given through 

a 

H NS ' a (v, = ± H q - a (-v, £) , 

H G (w,0=-ff 9 (w,0±^ < '(-«.0. 



(7) 



where the upper (lower) signs corresponds to the unpo- 
larized (polarized) case. The unpolarized singlet and the 
gluon have additional resonance-like contributions called 
the D-term [24,25] which ensures the correct symmetry 
and polynomiality properties (cf. eq.(4)). Note that the 
symmetries which hold for the matrix elements change 
for the H l s, due to the P ■ z, z factors in eq.(4). The 
unpolarized quark singlet is antisymmetric about v = 0, 
whereas the unpolarized quark non-singlet and the gluon 
are symmetric. The opposite symmetries hold for the 
polarized distributions. The hclicity flip GPDs, (labelled 
with Es), are found analogously (double integrals with 
respect to the Ks). 

Henceforth, we shall make the usual assumption that 
the ^-dependence of all of these functions factorizes into 
implicit form factors. One should bear in mind that in 
order to make predictions for physical amplitudes (for 
t ^ 0) these form factors must be specified. Note that 
the assumption of a factorized ^-dependence, as a general 
statement, must be justified within the kinematic regime 
concerned. It appears to be valid at small x and small t, 
from the HERA data on a variety of diffractive measure- 
ments. However, it appears to be wrong for moderate 
to large t and larger x [26]. The evolution used in this 
note assumes the same factorized i-dependence for each 
parton species, which therefore factorizes in the evolu- 
tion equations. Thus any effect of assuming a different 
i-dependence for quark singlet and gluon, which mix un- 
der evolution, is neglected (but could in principle be in- 
vestigated within the same framework). The effect on 
evolution of relaxing the assumption of factorization of 
the t-dependence altogether remains an open question. 

Taking the N — 1 moments (for even N) of the GPDs 
yields a sum of polynomials in the skewedness parameter 
£ [24]: 



dw N - 1 H(v,Z,t) 



N 
fe=0 



N 



1 dx f 1 F(x,y,t)x N - k y k . (8) 

■1 J-l+\x\ 



This powerful polynomiality condition arises because in 
indices of a local operator must be carried by 

either P or r (taking moments is equivalent to going 
from non-local to local operators). In building models for 
GPDs, eq.(8) always has to be fulfilled and thus strongly 
constrains any model for the GPDs. We incorporate it 
in a natural way in building our input distributions in 
section III. 



For the purposes of comparing to experiment it is nat- 
ural to define GPDs in terms of momentum fractions, 
X E [0, 1], of the incoming proton momentum, p, carried 
by the outgoing parton. To this end we adapt the nota- 
tion and definitions of [22] introducing two non-diagonal 
parton distribution functions (NDPDFs), T q and T q , for 
flavour, a: 



F q ' a [ X 1 
T"-" ( X 2 



(l-C/2) ' 
(l-C/2) 



(9) 
(10) 



where Vl G [-£,1],«2 G [-1,0 (see fig.(4) of [22]), 
£ = r + /p + is the skewedness defined on the domain 
C G [0, 1] such that £ = C/(2 - C) and C = x bj for DVCS, 
up to terms of 0(xbjt/Q 2 ). The inverse transformations 
between the vs and Xs is: 



X 1 - C/2 

Vl = -> T77T > v 2 



C/2 - X 2 



l-C/2 ' l-C/2 ' 

For the gluon one may use either transformation, e.g 



(11) 



F g (x,() 



(l-C/2)" 



(12) 



There are two distinct kinematic regions for the GPDs, 
with different physical interpretations. In the DGLAP 
[27] region, X > ( (\v\ > £), F q (X,C) and J*{X,Q 
are independent functions, corresponding to quark or 
anti-quark fields leaving the nucleon with momentum 
fraction X and returning with positive momentum frac- 
tion X — (. As such they correspond to a generaliza- 
tion of regular DGLAP PDFs (which have equal outgo- 
ing and returning fractions). In the ERBL [28] region, 
X < ( (\v\ < £), both quark and anti-quark carry posi- 
tive momentum fractions (X, ( — X) away from the nu- 
cleon in a meson-like configuration, and the GPDs behave 
like ERBL [28] distributional amplitudes characterising 
mesons. This implies that T q and T q are not indepen- 
dent in the ERBL region, indeed a symmetry is observed: 
T q (C,-X, C) = T q (X, () (which directly reflects the sym- 
metry of H q (v,(;) about v = 0). Similarly, the gluon 
distribution, T 9 , is DGLAP-like for X > ( and ERBL- 

like for X < (. This leads to unpolarized non-singlct, 
■pNS,a = T q.a_j:q.a^ and gluon qpq s wmcn are symmet- 
ric, and a singlet quark distribution T s = ^2 a j^i' a -\-^^ a 
which is antisymmetric, about the point X = C/2 in the 
ERBL region. Again the opposite symmetries hold for 
the polarized distributions. 



III. INPUT GPDS 

In this section we describe how to build input distri- 
butions, T q ' q ' 9 {XX,Qo), at the input scale, Qo, with 
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the correct symmetries and properties from conventional 
PDFs in the DGLAP region, for both the unpolarized 
and polarized cases. These input NDPDFs then serve as 
the boundary conditions for our numerical evolution. 

Factoring out the overall i-dependence we have for the 
quark for example (v\ >—£,): 



l-C/2 



f 1 r l ~\ x '\ F q - a (x' v') 

/ dx' dy >6( X >+Zy>- Vl ) V^Tlr- ( 13 ) 

J-l J-l+|x'| I- 1 



(l-C/2) 



Following [12,13,20] we employ a factorized ansatz 
which expresses the double distribution as a product 
of a profile function, it 1 , and a conventional PDF, /*, 

{i = q,g)- 



F q < a {x',y') = n q {x',y')f q < a {x') 
i {l-\x'\f-y' 2 

4 {i-\x'\y 

F°{x l ,y l ) = ^{x',y')f\x l ) 

15 ((l-M) 2 -y /2 ) 2 
16 (1-N'l) 5 

for quark of flavour a, where 



/V), (14) 



f(x) = xg(x, Qo)O(x) + \x\g(\x\,Q )e(-x) , 
p- a {x) = q a (x, QoMx) - (q a )(\x\ 1 Q )e(-x) . (15) 

The profile functions are chosen to guarantee the cor- 
rect symmetry properties in the ERBL region and their 
normalization is specified by demanding that the conven- 
tional distributions are reproduced in the forward limit: 
e.g. F 9 {X, ( -> 0) -> f 9 (X). We have explicitly checked 
this forward limit for our input codes against indepen- 
dent codes for the conventional PDFs (see fig.(l) later). 

So far, models for polarized distributions have been 
strongly based, for reasons of consistency, on closely as- 
sociated models for unpolarized distributions. For this 
reason we choose as our inputs these closely related pairs 
of models. Unfortunately, the polarized set of distribu- 
tions of Gehrmann and Stirling ('GS(A)') [15] is rather 
old now and is based on the MRSA unpolarized set [14] . 
Our other set is the very recent GRSV set from 2000 [17] 
which is based on the 1998 unpolarized model by GRV 
[16] (we will denote these GRSV00 and GRV98 respec- 
tively [30]). 

Having defined this model for the double distribution 
one may then perform the j/'-integration in eq.(13) using 
the delta function. This modifies the limits on the x' 
integration according to the region concerned: for the 
DGLAP region X > £ one has: 



F q > a {X,C,) = 



I 



+5 



dx'7T q 



Vl 



q a (x'). 



(16) 



For the anti-quark, since vi — —v\ one may use 
eqs.(10,13) with Vi,v\ — > — v\, and exploiting the fact 
that f q (x) = —q(\x\) for x < 0, one arrives at 



T q > a (X,() = 



dx TT q [ X , 



5°(|z'|). 



(17) 



The non-singlet (valence) and singlet quark combinations 
are given by: 

_ [H q < a (v l7 + H q ' a (-v u Q} 
l-C/2 

a 

= 2. T^cf2 ' (18) 

In the ERBL region (X > (,\v\ < £) integration over 
y' leads to: 

^' a (A,C)-|x 



I 



dx TT q \ X 



[ dx'ir q ( x 



T q - a {X,C) = — c * 



, Vl - X 



C 



q a (x')~- 
<f(\x'\) 



dx'lT q ( X ', 

dx'-K q I x 



q a (x') - 



I dx'n q (x',^f^-) 



<i a (W\) 



(19) 



This gives singlet and non-singlet distributions in the 
ERBL region containing four integrals over x' . The un- 
polarized singlet also includes the D-term on the right 
hand side of eq.(18) (in principle there is also an analo- 
gous term in the unpolarized gluon (D G in eq.(4)), but 
we choose to set this to zero, since nothing is known for 
the gluon D-term except its symmetry). We adapt the 
model introduced in [25] for the unpolarized singlet D- 
tcrm based on the chiral-soliton model: 



^ D (x,c) = e(c-x)^^-i,t = o) , 



(20) 



with D given by a truncated expansion in terms of odd 
Gegenbauer Polynomials: 



D(a) = (1 -a 2 )x 
-4.0C? /2 (a) - 1.2C 3 3/2 (a) - 0.4C 3/2 (a) 



(21) 
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This D-term is antisymmetric in its argument, i. e. about 
the point X = £/2 (in keeping with the anti-symmetry 
of T s , and H s about v = 0). It vanishes entirely in the 
forward limit and only assumes numerical significance for 
large £. This completes the definition of the input singlet 
and non-singlet quark distributions in both regions. The 
input gluon is defined along the same lines as the non- 
singlet. 



IV. NLO EVOLUTION OF GPDS 

The input GPDs, defined above, are continuous func- 
tions which span the DGLAP and ERBL regions, and 
evolve in scale appropriately according to the DGLAP or 
ERBL evolution equations, given below. Note that the 
evolution in the ERBL region depends on the DGLAP 
region (see eq. (25)) whereas the DGLAP evolution is in- 
dependent of the ERBL region. One remarkable feature 
is that the evolved functions remain continuous under 
evolution and continue to satisfy all appropriate symme- 
tries in the ERBL region. As can be seen, our numerical 
solutions manifestly exhibit these properties (see Section 
V), thus making us confident in the correctness of our 
numerical solutions. 

In the DGLAP region the singlet and gluon distribu- 
tions mix under evolution: 



dT s (y,(,Q 2 ) _ f 1 dz u (y A ^ s 



d\n{Q 2 



+ n-0/ *°<«.c,<f), 



with generalized DGLAP kernels [11]. The NS combina- 
tions do not mix under evolution so we omit them for 
brevity. 

The +-distribution for the DGLAP kernel, regulates 
the divergence as z — > y. In general it is defined as fol- 
lows: 

P(x,0+ =P(x,0~ S(x-1) f dx'P(x',0 (23) 

Jo 



plied to the integrals in eq.(22) is as follows: 



dz 



r(y,0[J s dzp( 



v y 



, y f c/y a \r ( y y 
+ cJo dzV { z cc 

In the ERBL region we have 



d\n{Q 2 ) 

y dz. 



k c 







c 



-V qq 



+ 1 



y z 
1 dz 



¥ 

C'C 



ldz vqq {u - 



c C 

e yq9 (H 

df G (y, C, Q 2 ) = r f 1 dz fy z 
dHQ 2 ) U v C VC'C 



+ 



f 



('-0 



dzV« [ V - - 



,hV "" 1 



(24) 



(25) 



with generalized ERBL kernels [11], where the bar nota- 
tion means for example, z/Q = 1 — z/C,. The upper signs 
correspond to the unpolarized case and the lower signs 
to the polarized case. 

The numerical implementation of the +- distribution, 
again applied to the whole kernel, takes the following 



The numerical implementation of the +- distribution ap- 
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form in the ERBL region: 




where the terms have been arranged in such away that all 
non-integrable divergences explicitly cancel in each term 
separately [31]. 



V. RESULTS 

We illustrate our polarized and unpolarized NLO 
GPDs in several ways, for two values of skewedness, 
C = 0.1, 0.0001, which are representative of the kinematic 
ranges accessible in HERMES and HERA kinematics, re- 
spectively. In particular, we quantify the effect of skewed 
evolution at NLO in the DGLAP region by comparing to 
conventional NLO DGLAP evolution of the same inputs, 
we present the relative influence of NLO skewed kernels 
by plotting the ratio of NLO to LO skewed evolution 
of common inputs and we also present the NLO GPDs 
themselves. In addition, we quantify the significance of 
the D-term and investigate the effect of changing the pro- 
file functions of eq.(14). 




MRSA' 
GRV 
GS 
GRSV 



- 0.9 - 

LL_ 

Q 
Q_ 

g 0-8 

CD 

0.7 f 

0.6 ' — — 

0.0001 0.001 0.01 0.1 

X 



FIG. 1. The ratio of the unpolarized and polarized input 
GPDs to conventional PDF inputs in the DGLAP region, for 
C = 0.0001. The forward limit GPD/PDF 1, for X > £ is 
observed in each case. The unpolarized GPDs are enhanced 
and the polarized GPDs suppressed close to X = £ = 0.0001. 

The skewed evolution of eqs. (22,25) of the input dis- 
tributions defined in Section. (Ill), for a particular £, is 
implemented by numerical integration on a fine grid in 
X and Q , using the CTEQ [32] evolution package, for 
initialization of the input and implementation of the Q 2 
evolution. The convolution piece was completely rewrit- 
ten. All code is written in Fortran77 and will be made 
available via the internet [23] . 

A. Comparison of skewed to conventional DGLAP 
evolution at NLO 




0.9 1 — — — 

0.0001 0.001 0.01 0.1 

X 



1 .05 




0.55 1 — — — — — 1 

0.0001 0.001 0.01 0.1 

X 

FIG. 2. Ratio of NLO-evolved GPD gluon over conven- 
tional PDF gluon at the input scale and several values of Q, 
for £ = 0.0001. The upper plot is for the unpolarized case 
(GRV98) and the lower one for the polarized case (GS(A)). 
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In fig.(l) we plot the ratio of our input gluon GPDs 
to the conventional input gluon PDFs on which they are 
based, in the DGLAP region (X > ( = 0.0001). The 
forward limit T{X > C, C, Qo)/f(X, Q ) -> 1, built into 
our input model of eqs.(13, 14, 15) by design, is illus- 
trated for each input PDF. In the region A « £, the un- 
polarized distributions MRSA and GRV98 exhibit the 
familiar enhancement (of about 15 — 20%) required by 
positivity constraints, T{X « (, (, Q )/f(X, Q ) > 1. In 
contrast, the polarized input model shows a suppression 
(of about 35 — 40%) in this region. We choose to plot re- 
sults for the gluon since this is most relevant to small- X 
exclusive diffractive processes at HERA, for which it has 
been argued that the inclusion of skewedness (at least 
at LO in the unpolarized case) leads to an enhancement 
(see e.g. [21,22,33]). 

In [21] it was demonstrated that LO skewed evolution 
in the unpolarized case yields a ratio consistently above 
unity with a strong enhancement near X = £ and a ratio 
tending to unity from above at large X ^> (. This is 
due to the probabilistic interpretation of the evolution 
at LO [34]. At NLO, this interpretation is lost in the 
unpolarized case [35] and thus the LO statement does not 
necessarily apply, since perturbative quantities become 
scheme dependent in NLO. 
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FIG. 3. Unpolarized NLO input and evolved singlet quark 
and gluon GPDs at small and large skewedness. The solid 
curves are the input MRSA GPDs at Qo = 2 GeV, the dotted 
ones show them evolved to Q = 10 GeV. The dashed curves 
are input GRV98 GPDs at Q = 2 GeV and the dashed-dotted 
ones show them evolved to Q = 10 GeV. The quark singlet 
is scaled by a factor of 10 -4 at £ = 0.0001 and by 10~ 2 at 
£ = 0.1. For £ = 0.1 the symmetry of the gluon GPD and the 
anti-symmetry of the singlet quark GPD are apparent about 
the point X = f/2 = 0.05. 

In the upper plot of fig- (2) we show the same ratio for 
the unpolarized evolved gluon distribution of GRV98 [36] . 
It is interesting and somewhat surprising to observe that 
this ratio drops progressively below unity even for X»( 
as Q 2 increases, and is as much as 6% below at X = 0.001 
for Q = 10 GeV. However, on inspecting the ratio of the 
NLO gg and gq kernels for the GPDs to the NLO gg and 
gq kernels of the forward PDFs, one observes that even 
for I > ( they can be very different from one another 
over a broad range of X, but will eventually settle at 
unity. An enhancement due to skewed evolution is still 
seen at NLO , but only very close to X — ( (for both the 
polarized and unpolarized cases). 

In the lower plot of fig- (2) we illustrate the polarized 
case using the GS(A) polarized gluon. Again we see the 
ratio drop progressively below unity for X > £. The 
size of the effect is similar, although this is difficult to 
see directly since the ratio of GPD to PDF is already 
smaller than unity at the input scale (for X ~ £) [37]. 



B. NLO GPDs and a comparison with LO 
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FIG. 4. Polarized NLO input and evolved singlet quark 
and gluon GPDs at small and large skewedness. The solid 
curves are the input GS(A) GPDs at Qo = 2 GeV, the dot- 
ted ones show them evolved to Q — 10 GeV. The dashed 
curves are input GRSV00 GPDs at Q = 2 GeV and the 
dashed-dotted ones show them evolved to Q = 10 GeV. The 
quark singlet is scaled by a factor of 10 -2 at £ = 0.0001. For 
f = 0.1 the symmetry of the polarized singlet quark GPD and 
the anti-symmetry of the polarized gluon GPD are apparent 
about the point X = </2 = 0.05. 

In figs. (3,4, 5) we present the input and evolved unpo- 
larized and polarized GPDs, respectively, for the singlet 
quark and gluon, and are encouraged to note that they 
are smooth, continuous functions of X. The lower plots 
in figs. (3, 4), with £ = 0.1, use a linear scale on the X axis, 
and hence illustrate very explicitly the symmetries in the 
ERBL region about the point X = £/2 = 0.05 which are 
manifestly preserved under evolution. These symmetries 
also hold for the upper curves (for £ = 0.0001) but are 
less visible on a logarithmic scale. These features indicate 
that our code is indeed correct and properly performs the 
NLO evolution in both the ERBL and DGLAP region. 
From fig. (3) and the upper plots of fig. (5) we note that 
both unpolarised input sets (GRV98 and MRSA') pro- 
duce similar results for the input and evolved GPDs. 

From fig. (4) and the lower left plot of fig. (5) we note 
with particular interest the result for the NLO evolution 
of the polarized quark singlet at small £. The GS(A) 
set assumes a rather different sea content to GRSV00, 
in particular the former has a positive quark singlet at 
the input scale and the latter a negative one. The NLO 
evolution keeps the GRSV00 distribution unchanged in 
shape, only altering its size, whereas the GS(A) distribu- 
tion is driven towards the same shape as well as overall 
sign as the GRSV00 distribution. This seems to indicate 
that the NLO generalised evolution favours the initial 
shape of the GRSV00 distribution and thus a negative 
polarized quark singlet at small £. At larger £, where the 
sea is of little numerical importance, the initial shapes 
are very similar and so one does not observe the same 
type of effect for GS(A). 




FIG. 5. Polarised ("A") and unpolarised ("V") singlet 
quark and gluon GPDs for ( = 0.0001 in the X > 10C re- 
gion. The input parton sets for each curve are the same as 
those given in Figs. 3 and 4. The quark singlet is scaled by 
a factor of 10 -4 in the unpolarised case and by 10 -2 in the 
polarised case. 

Although similar in shape (except for the polarized 
quark singlet at the input), the two sets of polarised dis- 
tributions differ considerably in size. The general evolu- 
tion effects which push partons to smaller X with increas- 
ing Q 2 , are naturally the same for both sets of inputs. 

In order to quantify the relative effects of including 
NLO skewed kernels in detail, we also evolved the same 
input GPDs using LO skewed kernels. In the upper and 
lower plots of fig. (6) we present the ratio of the evolved 
GPDs at NLO to those at LO, for the unpolarizcd and po- 
larized cases, respectively. Overall the effect of including 
NLO corrections is at most a 25% effect, within the kine- 
matic ranges illustrated. The regions furthest from unity 
in these ratio plots illustrate the ranges in X, relative 
to £, where the greatest effects are felt. The symmetry 
of the individual distributions about £/2 continue to be 
reflected in these ratio plots, which further indicates that 
the NLO evolution slightly alters the shape of the distri- 
butions relative to LO, but in such a way as to preserve 
the symmetries. 
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FIG. 6. The ratio of NLO to LO evolved GPDs for the 
unpolarized (upper plot, MRSA ) and polarized (lower plot, 
GRSV) cases. The solid and dashed-dotted curves shows the 
ratio at Q = 5 GeV and Q — 10 GeV , respectively. 

The spikes in the antisymmetric unpolarized singlet 
quark and polarized gluon ratios result from a slight loss 
of numerical accuracy in a few bins around the point Q/2 
(where the distributions have to pass through zero). The 
effect is magnified in this ratio. However, since the distri- 
butions vary by as much as twelve orders of magnitude 
(in the unpolarized quark singlet case), and the effect 
occurs close to a zero, this is not a cause for concern. 

The cusps observed at X = £ are due to the fact that 
the NLO over LO ratios approach the same value from 
the ERBL and DGLAP region very rapidly around the 
X = ( region. This change in shape in NLO over LO, 
which could be seen more clearly if the region around 
the X = ( were to be enlarged, is easily explained. In 
the LO kernels, rational functions of the type 1/(X — £) 



appear, which however combine to form a function which 
is finite at the point X = £. In NLO, functions of the type 
ln n (X -C)/(X- C)" 1 with n,nx = 0,1,2 appear, where 
again the rational functions cancel, but the logarithms do 
not. This leads to an integrable logarithmic singularity at 
NLO of a similar type as encountered in the NLO forward 
evolution, at the lower bound of the convolution integral 
in the RGE, y = Xbj, and leads to a strong enhancement 
of the NLO distributions in the region around X = ( 
as seen in the figures. The ratio of NLO to LO in the 
unpolarized case for GRV98 are very similar, within a 
few percent, to the plotted MRSA . In the polarized 
case the contrast is larger, until very large Q 2 , due to the 
markedly different input functions used for GRSV and 
GS(A) [38]. 



C. Influence of the D-term 



to 



FIG. 




x/C 



7. The influence of including the D-term on the 
GRV98-model quark singlet GPD at the input scale in 
the ERBL region X < for various values of £. The 
anti-symmetry which the D-term shares with the singlet, 
about X = £/2, is manifest by the symmetry about X — f/2 
in this ratio plot. 

Recall that in the definition of the unpolarized singlet- 
type quark input GPDs we included an extra piece in the 
ERBL region only, known as the D-term, which we mod- 
elled following [25] in eq.(21). This extra piece simply 
forms part of the definition of the input which is then 
evolved. In order to quantify its influence we consider 
in fig. (7) the ratio of the quark singlet with the D-term 
included to that with it omitted for various values of £. 
One observes that its numerical influence becomes pro- 
gressively less significant as £ decreases (at ( = 0.01 it is 
less than 1%). 
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D. Influence of the profile functions 

In section III we specified in eq.(14) specific cases of 
Radyushkin model [12,13] for GPDs based on double dis- 
tributions. The original model specified a more general 
form for the profile functions: 



216 



n(x,y) = 



r(26 + 2) [(1 - |a:|) 2 - t/ 2 ] 

2 2b + 1 r 2 (6 + i) (i - i^i) 2 ^ 1 



(27) 



So that eq.(14) corresponds to b = 1 for the quark and 
b = 2 for the gluon. It was suggested in [25] that vari- 
able b could be used as a fit parameter to extract GPDs 
from DVCS observables. It was also noted that the limit 
of very large b corresponds to the forward limit of £- 
independent GPDs, for all X > (. 

We have investigated this issue numerically and find 
that in the DGLAP region the GPDs tend to the PDFs 
very slowly as a function of b, particularly close to X = (. 
Meanwhile the behaviour of the GPDs in the ERBL re- 
gion is seen to change dramatically. We illustrate this 
in fig. 8 at small and large £ (i.e. ( = 0.0001,0.1, re- 
spectively), comparing b = 100 with the canonical val- 
ues for both the quark singlet and gluon, using GRV98 
forward distributions. As can be seen, in the DGLAP 
region (X > Q using such a large value of b has only a 
marginal effect. In the ERBL region a rather dramatic 
effect is seen at the input scale, particularly close to the 
symmetric point X = (/2. It is satisfying to observe 
(from the Q = 10 GeV curves) that under evolution both 
choices appear to tend to similar evolved distributions (as 
the memory of the details of the input distribution are 
washed out by the skewed evolution). We may conclude 
from the figures that if parameter b is to be used as a 
fit parameter to constrain GPDs using data, it would re- 
quire high statistics data on DVCS observables which are 
very sensitive to the details of the real part of DVCS am- 
plitudes in the ERBL region, such as the azimuthal angle 
asymmetry or the charge asymmetry (see e.g. [29]). We 
may also conclude that taking the very large b limit of 
Radyushkin's ansatz is not a particularly good way of nu- 
merically probing models in which the GPD is the same 
as the PDF in the DGLAP region, since the convergence 
appears to be very slow. 
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FIG. 8. Sensitivity of the input (Q = 2 GeV) and evolved 
(Q — 10 GeV) GPDs to the choice of the parameter b for the 
profile functions in eq.(27). We compare large b = 100 with 
the canonical values of eq.(14) (6 = 1.0 for the quark singlet 
and 6=2 for the gluon) , using GRV98 forward distributions 
at both small ( = 0.0001 and large £ = 0.1 skewedness (right 
and left hand plots respectively). The 6 = 100 case is shown 
with thickened lines. 



VI. CONCLUSIONS 

In this paper we have presented a complete numerical 
solution for the next-to-leading order (NLO) evolution 
equations for polarized and unpolarized generalised par- 
ton distributions (GPDs). We demonstrated that our 
solutions stay smooth and preserve all the symmetries 
of GPDs, giving us confidence in the correctness of our 
numerical implementation. We presented the relative ef- 
fect of moving from leading-order (LO) to NLO for a 
common input. We illustrate the effect of NLO skewed 
evolution on the magnitude and shape of two correctly- 
symmetrized input models based on conventional parton 
density functions, for two values of the skewedness pa- 
rameter, ( = 0.0001, 0.1, typical of HERA and HERMES 
kinematics. We demonstrate that in NLO, in contrast to 
LO, for small skewedness, the ratio of GPD gluon to PDF 
gluon in the unpolarized case does not have to be larger 
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than unity in the MS scheme employed here. This new 
effect is directly relevant to the study of small- a; exclu- 
sive diffractive vector meson and real photon production 
at HERA. The next step is to use these NLO GPDs to 
make NLO predictions for measurable physical processes 
[39-41], such as DVCS, at HERA [7,8], HERMES [9] and 
Jefferson Lab [10] . A detailed and direct comparison with 
all available data from a wide range of experimental ob- 
servables is ultimately required to establish accurately 
the details of the GPDs. 
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